%% Lindemann ratio simulation 
% (Go through all possible A and sigma)
%L_dot=centers of holes from voronoi analysis;
%d=length of each pixel at the unit of nm
%L_im= image for fitting
L_im=mat2gray(L_im); %normalize the data
a=sqrt(size(L_im,1)*size(L_im,2)*d^2/size(L_dot,1)*2/(3*sqrt(3)))*sqrt(3); % average hole-hole distance
x=[];
diffmat=[];
Array=(1:20)*0.003;
for A=Array 
  diff=[];
  for k=20:50
    sigma_ini=k*0.01*a; %simulate Lindemann ratio (k*0.01) from 0.2 to 0.5
    L_mat=zeros(size(L_im,1),size(L_im,2));  % simulated image
    for i=1:size(L_mat,1)
        for j=1:size(L_mat,2)
            Dis=(L_dot-[j,i])*d;
            L_mat(i,j)=sum(exp(-(Dis(:,1).^2+Dis(:,2).^2)/(2*sigma_ini^2)));   %intensity from every hole on a single point
        end
    end
    
    L_mat=A*mat2gray(L_mat);
    x=[x;20*0.01];
    diff=[diff,sum(sum(abs(L_im-L_mat)))]; %difference between experiment data and simulated results
  end
  diffmat=[diffmat;diff];
end

figure;imagesc(x,Array,diffmat);colorbar
hold on;
for i=1:size(diffmat,1)
 [~,ii]=min(diffmat(i,:));
 plot(x(ii),Array(i),'ro');
end
[min_value, min_index] = min(diffmat(:));
[min_row, min_col] = ind2sub(size(diffmat), min_index); %find the smallest difference point
plot(x(min_col),Array(min_row),'wo');xlabel('\sigma/a');ylabel('A')
sigma_ini=x(min_col)*a;
L_mat=zeros(size(L_im,1),size(L_im,2));
for i=1:size(L_mat,1)
    for j=1:size(L_mat,2)
        Dis=(L_dot-[j,i])*d;
        L_mat(i,j)=sum(exp(-(Dis(:,1).^2+Dis(:,2).^2)/(2*sigma_ini^2)));    
    end
end 
L_mat=Array(min_row)*mat2gray(L_mat);
figure;imagesc([0:size(L_mat,2)]*d,[0:size(L_mat,1)]*d,L_mat);axis equal;axis([0 size(L_mat,2)*d 0 size(L_mat,1)*d]);caxis([0 0.035]);
figure;imagesc([0:size(L_im,2)]*d,[0:size(L_im,1)]*d,L_im);axis equal;axis([0 size(L_mat,2)*d 0 size(L_mat,1)*d]);caxis([0 0.035]);